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We study spreading dynamics of nematic liquid crystal droplets within the framework 
of the long-wave approximation. A fourth order nonlinear parabolic partial differential 
equation governing the free surface evolution is derived. The influence of elastic distor- 
tion energy and of imposed anchoring variations at the substrate are explored through 
linear stability analysis and scaling arguments, which yield useful insight and predictions 
for the behaviour of spreading droplets. This behaviour is captured by fully nonlinear 
time-dependent simulations of three dimensional droplets spreading in the presence of 
anchoring variations that model simple defects in the nematic orientation at the sub- 
strate. 



1. Introduction 



Thin liquid films have surprisingly wide application in our daily life. From industrial 
coating and painting processes to printing, many current technologies require an under- 
standing of fluid flows in which one spatial dimension (the film thickness) is significantly 
smaller than the others (typically, the lateral scales over which film thickness changes). 
Under such circumstances one can use systematic asymptotic methods based on a small 
parameter (the representative aspect ratio of the film, which characterises the size of the 
free surface gradients) to simplify the full Navier-Stokes governing equations. Expanding 
the dependent variables of interest, such as fluid velocity, pressure, etc., in terms of this 
small parameter, one can obtain a much more tractable system of reduced equations 
for the leading-order quantities. Despite its simplicity, this approach has been used and 
experimentally tested many times, and has been very successful in describing the real 
physics in a wide range of flows. 

While plenty of work has been done with Newtonian fluids, this kind of systematic 
asymptotic treatment of flowing thin complex fluids, in particular liquid crystals, is still 



in its infancy (see Munch et al. ( 2006 ), Blossey et al. ( 2006 ) and Myers ( 2005 ) for examples 
of work on non-Newtonian, but not liquid crystal, thin film flows). Liquid crystals are 
anisotropic liquids, which typically consist of rod-like molecules. In a nematic phase, 
the rod-like molecules have no positional order, but they self-align to reach long range 
directional order. Therefore, to have a complete description of a nematic liquid crystal 
(NLC) flow, one needs to consider not only the velocity field, but also the orientational 



director field. In experiments on spreading nematic droplets, Poulard & Cazabat (20051 



found that NLC droplets spreading on a horizontal substrate exhibit a surprisingly rich 
range of instabilities, in the regimes where Newtonian droplets would only spread stably 
(see also Delabre et al. (20091 and Manyuhina et al. ( 2010| )). 
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As regards asymptotic (long wavelength) modelling of such flows, Ben Amar & Cum- 
|mings| ( |2001[ ) derived a strongly elastic model to describe the surface evolution of strongly- 
anchored NLCs, work that was extended by Cummings (2004) to the weakly- anchored 
case; while Carou et al. (2007) studied the model for blade coating of NLC in two di- 
mensional space in the limit of weak elastic effects. An alternative approach based on 
energetic arguments was presented by Mechkov et al. (20091. However, these different 
approaches lead to different predictions for the stability of a thin film; a discrepancy that 



(2001 


) and 


Cummings 


(2004) 



the same stress balance at the free surface of the film as in standard Newtonian flow, 
balancing pressure with capillarity. When this condition is modified to also include an 
elastic stress, results consistent with the energetics-based approach are obtained. More- 
over, these consistent results indicate that in the case of strong anchoring conditions at 
both, the solid substrate and the free surface of the nematic film, such a film is never 
unstable ( Lin et aZ.|2013 ) . Note, that in the weakly elastic limit of Carou et al. (2007) the 
effect does not appear at leading order and no effect on stability is seen. To summarise: 
many questions remain to be addressed regarding the instability mechanisms in free sur- 
face nematic flows. No fully consistent "lubrication" model for a three dimensional (3D) 
situation that can account for the weak anchoring effects that are crucial for instability 
has yet been proposed or studied. 

In this paper, we implement the long wave approximation to derive a model describing 
the three dimensional free surface evolution of a thin film of NLC on a rigid substrate. 
The model incorporates a novel weak anchoring surface energy formulation, and shows 
satisfactory behaviour in the vicinity of a contact line. Simple linear stability analysis 
permits mechanistic insight into how the anchoring energy influences the stability of a 
spreading NLC droplet. 



2. Model Derivation 

The main dependent variables governing the dynamics of a liquid crystal in the ne- 
matic phase are the velocity field v = (u,v,w), and the director field n — (^1,^2,^3), 
the unit vector describing the orientation of the anisotropic axis in the liquid crystal 
(an idealised representation of the local preferred average direction of the rodlike liquid 
crystal molecules). The director orientation is a function of space and time which, in 
the limit that director relaxation is fast relative to the flow timescale (the limit consid- 
ered here) is determined by minimising a suitably-defined total energy. Molecules like to 
align locally, a preference that is modelled by a bulk elastic (Frank) energy W, which is 
minimised subject to boundary conditions. In general, a bounding surface is associated 
with a given preferred direction for n; this preference is known as surface anchoring, 
and is modelled by an appropriate choice of surface energy. Anchoring can be tuned by 
appropriate treatment of a surface and may be either weak or strong. The stress tensor 
for the NLC is a function of the director orientation, hence elastic effects can strongly 
influence the fluid flow, giving rise to behaviour that differs markedly from the isotropic 
Newtonian case. 



2.1. Leslie-Ericksen Equations 

The flow of nematic liquid crystal may be described by the Leslie-Ericksen equations ( |Leslie 
1979). Neglecting inertia, and using over-bars to denote dimensional variables (dimen- 
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sionless variables will be without bars), the flow is governed by 
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representing energy, momentum and mass conservation, respectively. Here, A is a La- 
grange multiplier ensuring that the director n is a unit vector. The quantities W , G and 
II are defined by 

2W = K ((V • nf + |V A n\ 2 ) ; (2.4) 



G, 



-71 Ni -72 eikUk, e< 



1 



+ 



2 V dxj dx 
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UJ, 



dvj 
dx.j 



dx. 



n=p+w + ip„, 



(2.6) 
(2.7) 

where K is an elastic constant (this form of W (2.111 exploits the widely-used one- 
constant approximation ( De Gennes &: Prost|l995 )), 71 and 72 are constant viscosities; 
an over-dot denotes a material (total) time derivative; p is the pressure and ^ g is the 
gravitational potential. Finally, Uj is the viscous stress tensor, given by 

iij = ainkUpEkpniUj + a 2 Ninj + a 3 Njni + a 4 eij + a^e^nkUj + a 6 e~jkn k ni, (2.8) 

"3-^2, 72 = a 6 -a 5l 



where are constant viscosities (related to 7^ in Eq. (2.5) by 71 
and to each other by the Onsager relation, 02 + 0:3 = — a§). 

2.2. Nondimensionalization 
We make the usual long-wave scalings to nondimensionalize the governing equations 

(x,y,z) — (Lx, Ly,5Lz), (u,v,w) — (Uu,Uv,SUw), 



t 



L 



t, p 



p, W 



K 



W, 



(2.9) 



U'' r S 2 L^' S 2 L 2 
where L is the lengthscale of typical variations in the x and y directions, U is the typical 
flow speed; S = h /L -C 1 is the aspect ratio of typical variations of the film height h (a 
small slope assumption), and p — on/ 2 was chosen as the representative viscosity scaling 
in the pressure, since this corresponds to the usual viscosity in the isotropic case. Our 
choices of L and U are discussed later in the text. 

2.3. Energetics of director field 

It has been shown ( Ben Amar fc Cummings|2001| Cummings|2004 ) that with the above 
scalings, and p rovi ded the inverse Ericksen number K/{fiUL) = 0(1), the coupling terms 
in Eqs. (2.1 1- (2.3 1 between the energy and momentum equations, represented by G, 



can be neglected. The energy equations then reduce to the appropriate Euler-Lagrange 
equations for minimising the free energy of the film subject to the constraint n n = 1, 
corresponding to the limit of instantaneous relaxation of the director field. Imposing the 
constraint n ■ n = 1 directly we have a director field that is a vector on the unit sphere 
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characterised by two angles, 

n = (sin 9 cos (/>, sin 9 sin 0, cos 9), (2-10) 

for some functions 9(x, y, z, t) and (j>(x, y, z, t), which are the usual spherical polar angles. 
The leading order bulk elastic energy, under the long-wave scaling, is given by 

2W = 9 2 z + <j) 2 z sin 2 9 + 0(d). (2.11) 

The surface energy at the free surface z = h(x,y,t) is denoted by Q = Q(9) where 9 is 
the conical director orientation at the free surface, 

§ = e(x,y,h,t). (2.12) 

The surface energy Q takes its minimum when the director takes the preferred orien- 
tation 9 = 0. Within the long-wave approximation this corresponds to a director field 
perpendicular to the free surface: homeotropic surface anchoring. At the substrate z — 
we assume strong planar anchoring, 9(x, y, 0, t) = tt/2, with <f> specified. These anchoring 



assumptions are consistent with the experiments of Poulard & Cazabat ( 2005 ) (but not 
to all experimental spreading scenarios; in particular our model is not applicable to the 
case of fully degenerate planar anchoring at the lower substrate) . 

We carry out the free energy minimisation directly using a variational principle. The 
total free energy, J, consists of bulk and surface contributions. We write 

J= f f NWdSdz + f QdS, (2.13) 
Jo Jn Jn 

where fi is the domain occupied by the liquid crystal sample in the x-y plane and 
TV = K/(fj,UL) is the inverse Ericksen number. We consider the variations induced in 
J by small variations in the fields 9 and <j). The first variations must both vanish at 
an extremum and the sign of the second variations tells us whether or not we have an 
energy minimum. After an integration by parts, the vanishing of the bulk terms in the 
first variations of J leads to 

9 ZZ = ^ sin2(9 in Q U {0 < z < h}, (2.14) 

(cf) z sm 2 9) z = in ftU{0 < z< h}. (2.15) 

At the free surface, the surface energy Q is independent of the azimuthal angle <p (conical 
anchoring), hence a natural boundary condition on <f> emerges from the surface contri- 
bution to the first variation of J with respect to <p: <f> z sin 2 9 — on z = h. The angle cf> 
is thus independent of z; and with our assumption of strong anchoring at the substrate, 
we then have 

<f> = <f>(x,y) (2.16) 



determined by the imposed substrate anchoring pattern. For 9, Eq. (2.14) reduces to 



$zz = 0, and the strong planar anchoring condition is imposed on z — 0. We then have 

9 = a(x, yi t)z+^, (2.17) 

where a is determined by the condition that the surface contribution in the first variation 
vanish, 

Qs+Na = 0. (2.18) 
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2.3.1. Surface Energy 

For relatively thick films, the director angle 9 can easily adjust to the preferred values at 
each surface. As the film gets thin, and in particular near precursor layers or contact lines, 
there is a very large energy penalty to pay for bending between two fixed angles across 
a very short distance h. In this paper we assume the existence of a thin precursor film 
ahead of a bulk droplet, of thickness < b -C 1 (this is also the case in the experiments 
of Poulard & Cazabat (20051). To avoid a near-singularity in the director orientation 
within the precursor, we allow the anchoring to be relaxed as h — > b. 

To capture these two limiting behaviours for thick and very thin films, we propose that 
the change in director angle across the fluid layer, ah, approaches a prescribed value 8 
(the difference in the preferred angles at the free surface and solid substrate) as h — > oo; 
and approaches zero as the film thickness h —> b. Similar to the approach of Cummings 
\et al. (2011 ), we introduce an ad hoc anchoring condition based on specifying this change 



in director angle by 

ah = em(h), (2.19) 

where m{h) is a monotone increasing function of h with m(b) = and m(oo) = 1. With 
our assumption of homeotropic alignment at the free surface (and with the assumed 
lubrication scalings), = — tt/2. 

Though we did not specify the surface energy Q in the above, it is implicitly imposed, 
and easily recovered. Based on the above, the director angle 8 is g iven by 9 = (tt/2)(1 — 
zm(h)/h), so that the angle 9 at the free surface, defined by (2.12), is given as a function 
of h by 

§=*(l-m(h)). (2.20) 



The surface energy must satisfy equation (2.18). Since Eq. (2.20) is not trivially inverted 
to give h(9), we use the chain rule to obtain 



dG _ d9 _ m(h)m'(h) 
dh ~ d dh ~ h : 



(2.21) 



where Af — Q 2 Af. Equation (2.21 ) defines the surface energy g in terms of the film height 



h. The expression in terms of director angle at the free surface, 9, may be recovered by 



use of Eq. (2.20) 



Note that Eqs. (2.17) and (2.19) also imply that the bulk elastic energy in Eq. (2.11) 
becomes 

~2~~h?~ 



W = 



(2.22) 



As a result we have the contribution of nematic elasticity to the free energy as 



J = 





'Af 


2 "1 

m 


1 

la. 


2 


~h +g _ 



dS. 



(2.23) 



2.4. Momentum Equation 



For the momentum equations, Eqs. (2.2), balancing dominant terms gives 

OIL 

dx 



dtts 



dz 



an 

dy 



dt 2 3 
dz 



in dimensionless form. Based on the long-wave scalings, to leading order we have 



ti3 = (Ai+A 2 cos2(f))u z + A2Sm2(j)v z , t 23 = A 2 sin 2cf) u z + ( A x - A 2 cos 2<j>) v z , (2.24) 
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where A\ = ± + (0:5— a 2 ) cos 2 9+ct\ sin 2 9 cos 2 9+(a a +a 6 ) sin 2 9/2, A 2 = ol\ sin 2 0cos 2 9+ 
(a 3 + a 6 ) sin 2 9/2, and the cti are normalised by \i — a 4 /2. As a result, the leading order 
equations are 



-M9 Z 



Op 
dx 

9,2 



5_ 

c9.r 



{(Ai + A 2 cos 2(f>) u z + A 2 sin 2(f) v z } . 



{A 2 sin 2</> u z + (Ai — A 2 cos 20) i> z } , 



(2.25) 
(2.26) 
(2.27) 



where B — S 3 pgL 2 / \x\J is the Bond number. 

We assume that the normal component of the stress at the free surface balances surface 
tension (the isotropic component of the surface energy, 7) times curvature, and that the 
in-planc component of the stress is balanced by surface tension (surface energy) gradients 
in the plane of the surface. This yields the leading order boundary conditions: 

p + M9 2 z = -CV 2 /i, (2.28) 
-N{9 X 9 Z + 9 2 z h x ) + (Ax + A 2 cos 2(f) u z + A 2 sin 2(j> v z = NQ X , (2.29) 



N(9 y 9 z +9 2 z h y ) + A 2 sin2(j)u z + (A x -A 2 cos20) v z =NQ y , 



(2.30) 



where C = S 3r y/p,U is an inverse capillary number. Furthermore, using Eqs. (2.18)-(2.20), 
the equations of the tangential stress balances, Eqs. ( 2.29 1-( 2.30), reduce to u z = and 
v z =0. 



We solve Eqs. p^7| |- p^8| | for p 

p = B(h- z)-Na 2 -CV 2 h, 



(2.31) 



and substitute in Eqs. ( 2.25 )-( 2.26) to obtain u z and v z using the boundary conditions 
derived above: 



Du z = 

DVr.= 



(A\ — A 2 cos 2(f)) (p x + Naa x ) — A 2 sin 2(f) (p y + Naa y ) (z — h), 
(A\ + A 2 cos 2(f)) (p y + Afacby) — A 2 sin 2(f) (p x + Afaa x ) (z — h), 



where D — A 2 — A\. Finally, using conservation of mass together with the relations 



u dz 



u z (h — z) dz, 



v dz 



v z (h — z) dz, 



we obtain a partial differential equation governing the evolution of the film height: 



fi' + fa 



cos 2(f) sin 2(f> 
sin 2(f) — cos 2(f> 



V CV 2 h-Bh- 



= 0, 



where / is the identity matrix and 



fi = 



A x 



A\ A\ 



(h - z) 2 dz, j 2 = 



-Ao 



A\-A\ 



(h — z) 2 dz. 



(2.32) 



(2.33) 



Equations ( 2.32 )-( 2.33 ) represent a formidable analytical challenge. We simplify by 



approximating the integral expressions using the two-point trapezoidal rule, as 

2 + a 3 + a e a 3 + a 6 



fi 



h 



vh 6 



A 



(2.34) 



4(l + a 3 + a 6 ) 4(1 + a 3 + a 6 ) 

For — 1 < a>3 + a 6 < (which is the case for all common nematic liquid crystals), we have 
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A > v > 0. By including these quantities and our chosen surface energy Q from § |2.3.1[ 
the equation can be rewritten as 



where TV = 2 K/^UL and 

V = I AJ + /' 



h 3 V (CV 2 /i - Bh) + Af (mm'h - m 2 ) Wh 



= 0, 



cos 2(f) sin 2(f) 
sin 2cf) ~ cos 2(f) 



• V. 



(2.35) 



(2.36) 



2.5. Model Summary 



Our final model consists of the partial differential equation, Eq. ( 2.35 1, where V is defined 



in Eq. (2.36) with the anchoring condition at the free surface, m(h), and the anchoring 



pattern at the substrate, 4>{x,y), to be specified. We have five dimensionless positive 
parameters: A, v, C, B, Af, giving a solution space that is potentially very large. In the 
following analysis and simulations, we assume a balance between surface tension and 
gravity, setting C = B = 1, meaning physically that the typical length scale L considered 
is the capillary length, \/j/pg. Clearly, alternative choices can be made if a different 
balance, or different lengthscales, are to be considered. 



On the other hand, it has been shown (Lin et al. 20131 that the evolution equation 



for a NLC film in the limit of strong anchoring can be written in a variational or gra- 
dient dynamics form in line with such formulations for films of simple liquids (see, e.g., 
Mitlin| (|1993[); |Thiele| (|2010j )), films of mixtures ( |Thiele|[20lT| and surfactant-covered 
films ( Thiele et aL|2012 ) . We would like to point out that this is also true for the current 
model. In such a formulation, the evolution of the film thickness h follows a dissipative 
gradient dynamics governed by the equation 



Q(h)V 



SF 
~6h 



(2.37) 



where Q(h) is the mobility function and F is the free energy functional written as 

{Vh) 2 \ B 



F[h] = f 



C I 



-Ir 



dS + J. 



(2.38) 



(The contribution of nematic elasticity on the free energy functional, J, is given in 
Eq. ( 2.23[ ).) By introducing F into Eq. (2.37) and noting that the surface energy Q 
is coupled with the film thickness through Eq. (2.21), we obtain the evolution equation 



Q(/i)V [ -CV 2 h + Bh 



Afm 2 

YH 2 



(2.39) 



One should note that Eq. (2.35) and Eq. (2.39) are identical when Q(h) = h 3 



3. Analysis and Results 

In this Section we investigate some limiting cases of the model analytically, and carry 
out additional simulations for spreading nematic films and droplets in selected flow 
configurations. The time-dependent simulations that we report below are based on an 
Alternative-Direction-Implicit (ADI) method (as outlined by Witelski & Bowen (20031) 
with variable time stepping based on a Crank-Nicolson scheme; see Lin et al. (20126 1 for 
further details. 
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3.1. Influence of anchoring patterns at the substrate 

To gain some insight into our model we first compare two different uni-directional sub- 
strate anchoring patterns, in the simple case where flow is independent of y, and the 
fluid spreads uniformly in the x direction. Assuming <f) — (director orientation at the 
substrate parallel to the fluid flow direction) , Eq. ( |2.35 1 becomes 

h t + {\ + v) d x [h 3 {h xxx -h x )+N (mm'h - to 2 ) h x ] = 0. (3.1) 

If, on the other hand, we assume <j) = 7r/2, so that the substrate director orientation is 
perpendicular to the fluid flow (but still in the plane of the substrate), Eq. (2.35 ) becomes 

h t + (X- v) d x [h 3 {h xxx -h x )+N (mm'h - to 2 ) hJ = 0. (3.2) 



Equations (3.1 ) and A3. 21 are identical once time is rescaled by a constant. In this simple 



example, the substrate pattern only affects the flow timescale, effectively making the 
fluid viscosity (^-dependent. The NLC flows faster when the anchoring pattern is parallel 
to the flow direction (effective viscosity is smaller), and slower when it is perpendicular 
(effective viscosity is larger). 

To focus more on the influence of the anchoring patterns at the substrate, we next 
consider a weak conical free surface anchoring on 6 so that the director orientation is 
mainly determined by the strong planar anchoring at the substrate, i.e., 8 = ir/2 and 
m(h) = 0. The contribution of nematic bending elasticity then disappears. In particular, 



the integral expressions in Eq. (2.33) can now be evaluated exactly; there is no need to 
approximate them as was done to obtain Eq. ( 2.35[ ). The resulting equation is 



2h 3 ~ 

(V^h-h) 



= 0. 



(3.3) 



Figure [l] shows the solution of Eq. (3.3) computed using ADI-based simulations of a 
configuration where the anchoring imposed at the substrate appears as a striped pattern, 
as shown in Fig. [IJa): = ir/2 for x € (4n — l,4n + 1), n = 0, ±1 and ±2, and = 
otherwise. The initial film profile, shown as a surface contour plot in Fig. [ijb) , is taken 
as 

h(x, y, 0) = 0.45 tanh(-5(y - 10)) + 0.55, (3.4) 

with the front position being a straight line parallel to the x-axis, and spreading in the 
+y direction. The computational domain is defined by — L x x L x and y L y , 
with L x = 10 and L y = 20. The implemented boundary conditions are h x (±L x ,y,t) = 



h y (x,0,t) 



h xxx (±L x ,y,t) = h yyy (x,0,t) = 0, h(x,L y ,t) = b, where 



b is the thickness of the prewetting layer (precursor film) used to remove the contact 
line singularity. In the simulation scenarios that follow (in both this and the following 
Section) , the exact value given to b was found to have only a weak influence on spreading 
(with faster spreading for larger b) , but no influence of the exact value of b on the stability 
of the flow has been found. We use b = 0.1 here, and on uniform computational grids 
specified by the grid spacing 8x = 5y = 0.1; these values are found to be sufficient to 
guarantee numerical convergence. 

Figure [l|a) shows the evolution of the spreading fluid front. The anchoring pattern 
at the substrate clearly influences the spreading: as Fig. [lja) shows, in line with the 
observations discussed above, the front moves fastest when the fluid motion is parallel to 
the anchoring pattern, and slowest when flow is perpendicular to anchoring. As the an- 
choring changes periodically, the speed of the front transitions between the two extremes, 



giving rise to a sawtooth pattern. It should be noted that although Eqs. (3.1) and (3.2) 



are indicative of different wavespeeds for an isolated moving front, the amplitude of the 
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Figure 1. (Colour online) Spreading NLC film on a stripe-patterned substrate (A = 1, v = 0.5, 
M = 0). (a) The dashed (black) lines indicate the anchoring at the substrate. The solid (red) 
curves show the front position at At = 10 time intervals between successive curves. The initial 
front position is shown as a straight (red) line at y fa 10. (b - c) Surface contour plot of the film 
at t = 0, 50, respectively, (d) Cross section of the film, h(x = 0, y, t = 50). 

front perturbation does not increase linearly in time, as seen in Fig. [ija). Instead, the 
amplitude approaches a constant value. This constant amplitude is determined by the 
balance between different effective viscosities and surface tension effects. 

Figure [ljc) shows a surface contour plot of the profile at a late time, t — 50, and 
Figure [ijd) shows its cross section at x — 0. Note the absence of a capillary ridge 
behind the front, indicating the stability of the underlying flow. The sawtooth pattern 
that develops in the spreading front here, though reminiscent of a fingering instability, 
is no such thing: it is simply the result of the anchoring inhomogeneity imposed at the 
substrate. 

3.2. Influence of anchoring condition at the free surface 

We now analyse the effect of free surface anchoring on the director angle, 9. We begin 
by reviewing the linear stability analysis (LSA) of a simple flat film of height h = ho in 
the 2D case in which variations with respect to the y coordinate are neglected, so that 
both director field and flow are confined to the (a;, z)-plane ( Cummings et aL|2011 ). This 
analysis is found in practice to give a remarkably good indication regarding stability 
of spreading 3D droplets, considered in §3.3| below. In that section we present several 
simulations of stable/unstable 3D droplet evolution of a chosen surface anchoring function 
to illustrate the kind of behaviour that our model can reproduce in the presence of some 
simple substrate anchoring patterns. 

3.2.1. Linear Stability Analysis of a flat film 

With the director confined to the (x, z) plane, = 0, and no ^-dependence, Eq. ( 2.35 1 
reduces to Eq. (3.1). Assuming h = ho + £ and |£| <C ho in this equation, we find 

& + (A + v)h% K 



X 

m(h)m'{h)h — m(h) 2 



h 3 



(3.6) 



0, (3.5) 

where 

M(h) = 

By setting £ oc exp (ikx + ut), we obtain the dispersion relation 

u) = -(A + v)h 3 [k 4 + (1 - NM{ho))k 2 ] . (3.7) 

The flat film is thus unstable to sufficiently long- wavelength perturbations if MM {ho) > 
1. When this is the case, perturbations with wavenumbers k £ (0, k c ) are unstable, where 
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k c = ^AfM(ho) — 1 is the critical wavenumber. The fastest-growing wavenumber (for 
which the growth rate is the largest) is k m — k c j\j2, corresponding to the wavelength 

lm ~ ^(AfM(h )-l)J^' (3 ' 8) 

and the growth rate uj m — (A + v)hl{NM{h ) - l) 2 /4. 

3.2.2. Strong surface anchoring 

We consider firstly a strong homeotropic anchoring, given by m = 1. In this limit, the 
evolution equation becomes 

h 3 V {CV 2 h - Bh) - NVh\ = 0, (3.9) 

Note that, the elastic contributions to the governing equation are purely diffusive. A 
version of this limit was derived via alternative energetic considerations by |Mechkov| 
et al. (2009) (see also Lin et al. ( |2013[ ) for a more in-depth discussion of the strong 



anchoring case). Although the director field corresponding to strong anchoring becomes 
singular as the film height h — >• 0, the PDE governing the film height in this limit is well- 
behaved and will never exhibit an instability. This observation suggests that the weak 
free surface anchoring, necessary on physical grounds for the director to be nonsingular 
as the film height goes to zero, is key for the instability mechanism. 

3.2.3. Weak surface anchoring 

For a weak surface anchoring, there are many possible forms for m(h) that satisfy our 
basic requirement m(b) = 0, m(oo) = 1. Here, as an example, we take 

where a and j3 are positive constants that tune the relaxation of the anchoring for 
film heights larger than the precursor, and f(h; b) provides the "cutoff" behavior as the 
precursor is approached. In the simulations presented in this paper we choose f(h; b) = 
[tanh((/i — 2b) /w) + l]/2, where w fixes the size of the /i-range over which m(h) is turned 
off as h — > b (w — > gives a simple discontinuous switch; we assign a small positive value, 
w = 0.05, to smooth this behaviour). This choice for m(h) ensures that the director 
field for thin films, roughly less than b, lies in the plane 9 = tt/2, with (j> dictated by 
the substrate anchoring conditions. We note that the exact functional form given to 
m(h) does not influence the results to any significant degree, as long as m(h) changes 
sufficiently rapidly for h ~ b. 

Figure [2] (a) shows the anchoring condition at the free surface, Eq. ( 3.10[ ), with 



a 



(3 = 1 and b = 0.1. It can be seen that the anchoring condition approaches for thin films 
and increasingly approaches 1 when the film thickness gets thicker. Figure [2] (b) shows 
the function Nm{h) — 1 for N = 0.2 (dashed (black) curve) and for N = 1 (solid (blue) 
curve). As demonstrated in j 3.2. 1| Eq.(3.7), the flat film with thickness ho is unstable 



if Afm(ho) — 1 > 0. One can see that for = 0.2, flat films are always stable while for 
Af = 1, there exists a range of film thicknesses (the critical values are marked by circles 
(red) in Fig. |2](b)) that exhibit instabilities. 



3.3. Numerical results 



In this section we present several simulations of three dimensional spreading nematic 
droplets, in which the influence of the anchoring condition at the substrate can be directly 
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(a) m(h) 




0.3 . 0.6 

(b) NM(h) - 1 



Figure 2. (Colour online) (a) The anchoring condition, m(h), as defined in Eq. (3.101 with 
a = 8 — 1 and b — 0.1, w — 0.05. (b) NM(h) — 1 as a function of h where M(h) is defined in 
Eq. (3.61. The solid (blue) curve is for Af = 1 and the dashed (black) curve for Af = 0.2. The 
red circles indicate the critical values where AfM(h) — 1 = 0. 



investigated. In particular, as test cases we consider spreading on substrate anchoring 
patterns that mimic the director structure near different types of nematic defects, in order 



to analyse the influence of local director structure on spreading (see Lin et al. (2012a) for 
a related analysis of strictly 2D "defects"). Such defects are classified according to their 
topological winding number s: as a small planar circuit around the defect is traversed 
exactly once, the director field rotates through an angle 27rs. Specifically, we have 



4>(x, y) = s tan 1 ^— ) 



(3.11) 



While this description cannot capture the true physics very close to the defect centre 
(the director field notion breaks down there and one has to introduce a tensorial order 
parameter for a detailed description; see e.g. De Gennes & Prost ( 1995 ) for a discussion) 



we suggest that it may provide a reasonable description of the macroscopic free surface 
evolution in the presence of a pinned defect. More importantly, for different choices of 
s, Eq.(3.11) provides examples of generic surface anchoring patterns that provide useful 



demonstrations of our model behaviour. 



We solve numerically Eq. ( 2.35 ) on such an anchoring pattern, our choice of parameters 



guided by the LSA results presented in § |3.2.1| above. The computational domain in all 
cases is chosen as — L x ^ x ^ L x and — L y ^ y ^ L y with L x = L y = 20, and with the 
boundary conditions 



h(x, ±L y , t) = h(±L x , y, t) = b, h y (x, ±L y , t) = h x (±L x , y, t) = 0. 



(3.12) 



For our first set of simulations we take as initial condition a smoothed cylinder of radius 
10 and height ho — 0.2, with a precursor film of thickness b = 0.05 covering the rest 
of the domain. While the exact functional form for h(x, y, 0) only weakly influences the 
consequent evolution, we give it here for defmiteness: 



h(x, y, 0) = tanh (-(v^+7 - 10)) + 



Finally, the parameters appearing in Eq. (2.36), (3.10) are chosen as a = 1, j3 



(3.13) 
T A 



1 

and v = 1/2. The value of Af is given in each figure caption. 

Figure [3^b) and (c) shows the evolution of a stably spreading nematic droplet, for 
Af = 0.2. The anchoring at the substrate mimics the director structure near a defect of 
type s = —1/2, shown in Figure ^ a). Due to the non-uniformity of this pattern in the 
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(a) anchoring condition (b) t = (c) t = 10000 

Figure 3. (Colour online) Spreading NLC droplet for A/* = 0.2 and s = —1/2. (a) The anchoring 
condition at the substrate, (f)(x,y), as defined in Eq. (3.111. (b) The initial condition at t = 0. 
(c) The droplet evolution at t = 10000. 
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Figure 4. (Colour online) Spreading NLC droplet for N = 1 and s — —1. The initial condition 
is the same as in Fig.|3|b). The anchoring condition at the substrate is shown in (a). The droplet 
evolution at t = 500 and t = 1000 are shown in (b) and (c), respectively. 



radial direction, the droplet spreads asymmetrically, but in the manner that might be 
expected for the prescribed anchoring pattern. Consistent with the results of our LSA 
in § 3.2.1 we do not observe any free surface instabilities: the corresponding flat film is 
stable for the chosen parameters. 

Figure [4] shows the evolution of a spreading nematic droplet, for J\f = 1. The anchoring 
at the substrate mimics the four-fold symmetric director structure near a defect of type 
s = —1, shown in Fig.Qa). Again, due to the non-uniformity of this pattern in the radial 
direction, the droplet spreads in the manner that might be expected for the prescribed 
anchoring pattern. In addition, we observe rich pattern formation on the droplet surface, 
as illustrated by Fig. [4] (b - c) . 

By the analysis of § |3.2.1 a flat film of the same thickness ho = 0.2 is unstable for 
Af = 1. The most unstable wavelength for this case is l m « 27T, predicting that there will 
be about 3 humps on the droplet surface for the chosen initial condition. We find that 
the results of simulations are consistent with these predictions, see, e.g., Fig. Hb). The 
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(a) anchoring condition (b) t = 5 (c) t = 30 

Figure 5. (Colour online) Spreading NLC droplet for AT = 10 and s — 1/2. The initial condition 
is the same as in Fig.|3|b). The anchoring condition at the substrate is shown in (a). The droplet 
evolution at t = 30 is shown in (b). 




cross section in the radial direction of this figure shows 3 hump-like structures, specified 
by (in 3D) one raised ring with one spherical hump at the centre. 

In order to confirm that this comparison between the LSA and nonlinear simulations 
extends to other parameter values, we next consider JV = 10. Here, the LSA predicts 
that the most unstable wavelength is /,„ « 1.6 and the number of humps for the drop 
considered should be about 12. Figure [5] shows this case, for a droplet spreading on an 
anchoring pattern given by Eq. ( 3.11[ ) with s — 1/2. We again find remarkably good 
agreement between the LSA prediction for the flat film, and the observed simulation 
for the cylindrical droplet. Figure [5jb), for example, shows that there are 6 rings that 
form, corresponding in the cross section to 12 humps. Also note that the time scale 
for instability development is much shorter for Af — 10. Here, the unstable pattern has 
developed already at t = 5, see Fig. gb), while for Af= 1 we have to wait until t — 500, 
as shown in Fig.^b). This finding is also consistent with the LSA predictions. We remark 
also that the type of structures seen in Fig. [5^c) are remi niscent of certain free surface 
structures seen in the experiments of Poulard & Cazabat ( 2005 ) (albeit within a much 
more complicated setting in those experiments). 

Finally, in Fig. [6] we show a spreading NLC droplet on a radially symmetric anchoring 
pattern, which mimics the director structure near a defect of type s — 1. The parameter 
Af is set to unity. Instead of using a cylindrical cap as the initial condition, here we choose 
one that is close to a spherical cap, shown as the solid (black) curve in Fig. |6jb), a con- 
figuration for which our linear stability analysis is certainly not applicable. Nonetheless 
we note that, taking h in the LSA to be the mean initial droplet height, one might 
anticipate instability for these parameter values. As the droplet spreads, the front re- 
mains circular for all time, while the surface exhibits radially symmetric instabilities, as 
anticipated. The instabilities appear as a ring (two humps in the cross section), which 
eventually closes up into a single central hump for long times. 



4. Conclusions 

We have presented a new model that describes three-dimensional spreading of thin 
films and droplets of nematic liquid crystal. To the best of our knowledge this is the 
first model of this kind to account for the effect of director variation in three dimensions 
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Figure 6. (Colour online) Spreading NLC droplet for N — 1 and s = 1. The anchoring condition 
at the substrate is shown in (a). The cross sections of the droplet, h(x,y = 0,t), are shown in 
(b) for t = (solid (black) curve), t = 300 (dotted (blue) curve) and t — 600 (dashed (red) 
curve). The contour of the droplet at t = 0, t = 300 and t — 600 are shown in (c), (d) and (e), 
respectively. 



on the shape of the overlying free surface. The stripe pattern of nematic films in a 3D 
setting was already analysed by Sparavigna et al. ( 1994 ), Lavrentovich & Pergamenshchik 



(1994), and Manyuhina & Ben Amar (2013) assuming the film remains a flat film. There 



the predicted instability mechanisms depend on the ratio of the various elastic constants 
while here the presented mechanism results from the coupling of free surface modulations 
and director orientation as described by the one-constant approximation. 

Strong anchoring boundary conditions on the director at both boundaries are not 
suitable to describe a very thin spreading film (the director polar angle 6 becomes singular 
at a contact line, leading to a strong diffusion, which is always stabilising). Instead, we 
impose weak conical surface anchoring on the polar angle 9, with the anchoring energy 
given by Eq. (2.21). The anchoring at the substrate z = is taken to be strong and 



planar, with the azimuthal director angle 4>(x, y, 0) specified. Our formulation preserves 
the property of strong anchoring when the film is thick, while allowing the director to 
relax to a state of planar alignment (though with anisotropy entering through nonuniform 
azimuthal patterning) when the film is very thin. The resulting equation for the film or 



droplet evolution is a fourth order nonlinear parabolic PDE, Eq. (2.35) 



A simple linear stability analysis of Eq. (2.35) in the case of purely 2D flow predicts 
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that a fiat film may be unstable under certain conditions. The strong anchoring limit leads 
to a purely diffusive contribution from the elastic effects that always acts stabilising; but 
weak anchoring can lead to instability. The physical mechanism is based on a coupling 
of the degrees of freedom of director orientation within the film and at its surfaces and 
the shape of the free surface itself. Even in the limit of instantaneous director relaxation 
considered here this coupling gives rise to an instability mechanism active in the film 
thickness range where anchoring and bulk elastic energies compete. A second mechanism 
that can lead to patterned spreading (which might be viewed as an instability in the 
advancing front) is that the anchoring condition on the azimuthal angle ((f)) at the solid 
substrate affects the speed of spreading. A drop spreads faster/slower when the substrate 
anchoring is parallel/perpendicular to the flow. For a substrate characterised by non- 
uniform anchoring conditions, the fluid front advances non-uniformly, in line with the 
prescribed anchoring patterns. This behaviour is exemplified by the analysis of a film 
spreading over a substrate with a striped anchoring pattern, leading to evolution with a 
sawtooth pattern in the advancing front, shown in Fig. [I] 

We carried out numerical simulations of 3D spreading droplets for a variety of substrate 
anchoring patterns, focussing particularly on patterns that mimic the director structure 
near topological defects. Our simulations (including more than are reproduced here) 
indicate that (i) the flat film stability analysis serves as a remarkably good indicator of the 
stability of more complex spreading droplets, provided that the initial "droplet height" is 
well characterised; and (ii) although substrate anchoring clearly affects spreading speed 
and the shape of the spreading front, it does not appear to influence the global free 
surface stability of spreading droplets. 

Though simplified, the proposed model and the reported simulations provide valuable 
insight into the dynamics of spreading nematic droplets and films as observed experimen- 

~ (|2005| ), |Delabre et aL| ( |2009| and |Manyuhina et al.\ ( f2010[ ). 



tally by |Poulard fc Cazabat 



The model as given by Eq.(2.35) is rather general, relying only on the validity of the 



lubrication scaling (which in turn relies only on the droplet aspect ratio); the strong 
anchoring condition at the substrate; and the two-point trapezium rule approximation 
for the integral expressions appearing in Eq. ( 2.32[ ). Note that we propose and use a 
particular reasonable form for the anchoring function m(h) only where it is necessary to 
carry out simulations or to demonstrate possible (in)stability regions. Thus, whenever 
the anchoring function m(h) is obtained experimentally (as an empirical function to be 
fitted) Eq. (2.35 ), is applicable, and its predictions for the stability of a suitably thin flat 
film should be valid. 

However, there is still much to be done in order to elicit the full story in all its com- 
plexity. The results presented here, in particular regarding the influence of substrate 
anchoring patterns, clearly represent only a small subset of the possible spreading be- 
haviour. Only very simple spreading scenarios and anchoring conditions are studied here, 
and it would clearly be of interest to simulate droplets spreading over more complex sub- 
strate patterning; for example, droplets spreading over several model defects as might be 
relevant in physical experiments. Our suggestion that the proposed substrate anchoring 
patterns may be thought of as idealised representations of defects in physical flows may 
of course also be questioned: it is known that the continuum nematic description used 
here breaks down in a small (nanometers) region around any defect, so our model cannot 
give an accurate description within such a defect core. Nonetheless, our simulations give 
some useful insight as to the effect that patterned planar anchoring can have on droplet 
evolution; and the similarity of Fig.[5|c) to parts of Fig. 2(c) in Poulard & Cazabat (20051 
is intriguing. 

While qualitatively illuminating, some aspects of our model are undoubtedly overly 
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simplistic: our contact line regularisation may not adequately model the true physics as 
the ultra-thin precursor film is approached and molecular effects such as van der Waals' 
interactions become important; and indeed it may well not capture the true behaviour of 
the anchoring conditions. The effect of finite surface anchoring energy at the substrate 
may also need to be taken into consideration: the relative anchoring strengths at two 
bounding surfaces were found to be important in the transition behaviour of the director 
field (albeit in the presence of an applied electric field) by Barbero & Berberi ( 1983 ) (see 
also citing works) . In future work we plan to introduce improved models for these aspects 
of the problem. 
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